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Static properties and Stark effect of the ground state of the HD molecular ion 
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We have calculated static properties of the ground state of the HD + ion and its lowest-lying 
P-state without making use of the Bom-Oppenheimer approximation, as was done in the case of H 2 + 
and D 2 + [Phys. Rev. A 58, 2787 (1998)]. The ion is treated as a three-body system whose ground 
state is spherically symmetric. The wavefunction is of generalized Hylleraas type, but it is necessary 
to include high powers of the intemuclear distance to localize the nuclear motion. We obtain good 
values of the energies of the ground S-state and lowest P-state and compare them with earlier 
calculations. Expectation values are obtained for various operators, the Fermi contact parameters, 
and the permanent quadrupole moment. The cusp conditions are also calculated.The polarizability 
was then calculated using second-order perturbation theory with intermediate P pseudostates. Since 
the nuclei in HD + are not of equal mass there is dipole coupling between the lowest two rotational 
states, which are almost degenerate. This situation is carefully analyzed, and the Stark shift is 
calculated variationally as a function of the applied electric field. 

I. INTRODUCTION 

By observing Rydberg states of HD, it should be possible to extract properties of the HD + ion, 
which is the core of the excited molecule. Such experiments [1] have been performed on Rydberg 


states of H 2 and D 2 giving highly accurate polarizabilities of H 2 and D 2 ions. Our recent calculations 
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[2] for the polarizabilities of the spherically symmetric ground state of the homonuclear molecular 
ions agreed very well with the high precision experimental results [1]. We showed that it is possible 
to avoid the Bom-Oppenheimer approximation (BO), which uses a product of electronic and nuclear 
wave functions, by following the same approach in diatomic molecular systems as has traditionally 
been used in two-electron atomic systems [3]. We simply treat the molecular ion as a three-body 
quantum system. 

We showed in our previous paper [2] that it is possible to modify the generalized Hylleraas 
basis functions, which had worked very well in a number of two-electron systems, to apply to two 
heavy particles and one light particle, thereby treating all three particles on an equal footing. Since 
no approximations are made (beyond the finite expansion length) in the calculation of eigenvalues 
by a variational principle, the results obtained, which include the kinematic effects of * rotation and 
“vibration” and electron excitation, should be better than in the BO approximation, which assumes 
that there is no direct coupling between the electronic and nuclear motions. These effects are 
automatically included when we calculate three-body wave functions. 

It is well known that the power-exponential form of a Hylleraas function is well suited to 
describe the motion of electrons, which are not localized. In order to describe the relative motion of 
heavy particles which are fairly tightly localized, it is necessary to modify the forms of the functions 
by introducing very high powers of the distance between the heavy particles. This significantly 
improves the convergence of the energy values. Using the wave functions obtained by minimizing 
the energy using the Rayleigh-Ritz variational principle, just as for H 2 + and D 2 + [3], we calculate 
expectation values for various operators, the Fermi contact parameters, the permanent quadrupole 
moment, and the cusp conditions. The energies obtained are compared with those given by Moss [4]. 
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II. FORMULATION 

The unperturbed Hamiltonian of the HD + system is 


H=-— V*-V*+2[ 


M T 


M r 




K- R d I 


( 1 ) 


where R p , R D and R e are the position vectors of the proton, deuteron and electron, respectively, 
measured in units of the Bohr radius a 0 , M P and M D are the masses in units of the electron mass, and 
the energy is in rydbergs. The dipole operator d is defined as 

d=2(R p +R D -R e )t, (2) 


where Sis a unit vector in the direction of the applied electric field. 

We make the usual transformation to the center of mass system: 


Rcm 


M p R p +M D R D +R e 

M p +M D +l 


^ 1,2 Rpj) Rf 


(3) 


In terms of these new coordinates, excluding the center-of mass contribution, the Hamiltonian takes 
the form 


H=-—V\ V 2 2 -2V.-V, + 

t* D 


r r r 2\ 


2_ 
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(4) 


and the dipole operator becomes 
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d = 2 


R CM + 


M t 


r . + - 


M p +1 




■e. 


(5) 


where 


M p 

M T -M p +M D + l, ^p~~^ + j ’ 


and n D = 


Mp 

AVI’ 


( 6 ) 


(Notice that in this case the mass polarization cross-term in // is of the same order as the other 
kinetic-energy operators, rather than being small as it is in a two-electron atom like helium.) 
Numerically, Mp= 1 836.152 701 and M D = 3 670.483 014. It is worth emphasizing again that we are 
treating all three particles on an equal footing and will generally not refer to any special “molecular” 
quantum numbers except perhaps to clarify the physical situation.. 


III. CHOOSING THE FUNCTIONS 

We are interested in calculating the static properties of the L=G ground state of HD + including 
the binding energy of the ground state, and in order to compute the dipole polarizability it is also 
necessary to construct a set of pseudostates of orbital angular momentum L=l. In the past [5] we 
have used the following types of Hylleraas functions successfully for these purposes: 




V* 


-ar,-br n -a r 


lmn=0 


l m n 
1 r 2 r 12 


(7) 
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and V p = -cos(0 t /2)(f+g)Dl * -sin(0 |2 /2)(/'-£)D 1 l (8) 


where the radial functions f and g are given by 


Kr v r 2 ,r n )=e 


% 

a P r iz r y' j r ‘J * 

r \ a ijk r 1 r 2 M2 


(9) 


*/Jt=0 


9 > 


and g(r 2 ,r 1 ,r 12 )=e' Y(r,+r2) ^' (X ^ 1 V 2 £ d ijk r^r{r^. 


ijk=0 


( 10 ) 


The rotational harmonics D, +1 and Df 1 used in Eq. (8) are those defined by Bhatia and Temkin [6], 
and r, 2 =r, -r 2 . In this case there is no exchange between the nuclei and therefore two independent 
functions are required in Eq. (8) because the electron can be in an L=1 state with respect to the 
proton or with respect to the deuteron. For convenience in optimization, the nonlinear parameters 
have been kept the same in the two functions f and g but they could have been chosen different in 
the two functions. 

Although the exponentials and powers appearing in Eqs. (7), (9), and (10) are efficient for 
describing the correlations between the nuclei and the electron, we know from the BO model (which 
is certainly a good first approximation) that the intemuclear motion will not be well described that 
way. It is much better described by Gaussian-like functions centered around the equilibrium 
positions of the nuclei. The problem is to adapt our trial functions so that they describe the motion 
of the nuclei satisfactorily. The method we use is to increase powers of r 12 appearing in Eqs. (7), (9), 
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and (10). This has been illustrated for H 2 + in Fig. 1 of Ref. [2] which shows that the ground state 
wave function obained from the Morse potential [7] can be well approximated by the form r i2 e >f '~ , 
provided N is large and a s =N/2. This leads us to make a simple modification of the trial functions 
of Eqs. (7), (9), and (10): multiply Eq. (7) by r 12 ° and Eqs. (9) and (10) by r l2 p . 


IV. PROPERTIES OF HD + 

The energy eigenvalues are calculated by the usual Rayleigh-Ritz variational principle which 
requires minimizing E with respect to variation of the nonlinear parameters in the expression 


CP|T> ' 


(ID 


Here Y can be Y 0 or Y p giving E 0 or Ep. For Q 0 = 13 (560 terms in the expansion), N 0 =12 we obtain 
Eq = -1.195 795 889 Ry which compares well with -1.195 795 931 Ry obtained by Moss [4], 
differing by only 4.2xl0' 8 Ry. Similarly, for L=l, for Q p =11 (728 terms in Eq. (9)), N p =l we 
obtain the lowest eigenvalue E p = - 1.195 371 602 Ry which is fairly close to -1.195 396 256 Ry 
obtained by Moss [4], but differing this time by 2.5xl0‘ 5 Ry. 

Using the wave functions obtained variationally, expectation values of various quantities are 
calculated for the ground state as well as the lowest L=1 state. They are given in Table I. In the table, 
r represents the distance of the electron from the center of mass of the nuclei and z is the projection 
of r onto r, 2 and they are given by the expressions 
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^ 2 _ (l+c 3 )r t 2 Hl-c 3 y 2 2 (I-C3V12 
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and Q=^-(r?2 + r 2 ~3z 2 ) is the permanent quadrupole moment. In Eqs. (12) and (13) c 3 = — . 

2 M p +M D 

In the Table energies are in units of Ry, and distances are in units of a & The cusp values with respect 
to r, and r 2 are close to /x P and fx D , respectively. The Dirac 6 functions 6(r,) and 6(r 2 ) required in 
the calculation of the Fermi contact term are also given in the Table. There does not seem to be 
another comparable calculation of these properties. 


V. POLARIZABILITY 

Having established that accurate wave functions can be obtained for the ground as well as 
the excited states, we can now calculate the dipole polarizability a, by ordinary second-order 
perturbation theory: 


vE (o| yy |o) v 


V £ 0 


(14) 


where d is the dipole operator defined in Eq. (5), and we can drop the center-of-mass coordinate 
which has no direct effect on the polarizability. For this purpose we use the ground state wave 
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function obtained with 560 terms and for intermediate states we keep N p = 7 as we vary Not 
surprisingly, we find that almost all of the calculated polarizability comes from the lowest-lying p- 
state in the intermediate sum in Eq. (14). This is due to the very small difference 
A 0) = 4.24287 x 10 ' 4 /?y between the ground and first excited states. Because the polarizability is so 
sensitively dependent on this difference we have chosen to replace our calculated value by the more 
accurate value given by Moss, A Q1 =3.996750x 10’ 4 Ry.The polarizability due to this single term is 
then 392.0814 a 0 3 , and the sum over all the remaining intermediate states contributes 3.2076 a 0 3 . The 
convergence of this second component (cT) and the energy of the lowest p-state are displayed in 
Table II. If the ground state is excluded, there is a corresponding modified spherical polarizability 
of the lowest p-state (a,,) due to coupling only with s-states, not d-states. The convergence of this 
quantity and the energy of the ground state is shown in Table HI. The total polarizability, 
a = 395.289 a 0 3 , is very much larger than might be expected from the results [l,2,][8]in 
H 2 * and D 2 .This can be traced to the fact that HD + is heteronuclear, and so there is coupling 
between the ground state and the very low-lying first p-state, which would not happen in the 
homonuclear cases. We will discuss this situation in detail in the Appendix. 

It is essential to remember that the perturbation formula is valid only if the perturbing 
potential is small compared to the spacing of the unperturbed energy levels. In the homonuclear case 
this condition is easy to satisfy, since the lowest p-state is not coupled to the ground state, but for 
HD + where there is coupling to the low-lying p-state the calculated polarizability is meaningful only 
for extremely weak electric fields, much weaker than those found in Rydberg states of the HD 
molecule. In the next section we will quantify this statement and suggest a way to compute the Stark 


effect nonperturbatively. 
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VI. NONPERTURBATIVE CALCULATION OF THE STARK EFFECT 

To quantify the applied electric field in terms that are scaled conveniently for the present 
problem, let the perturbing potential be that of a fixed point (electronic) charge at a large distance 
R from the molecular ion, so that only the dipole term needs to be retained: 


V{R) = — , 

R 2 


(15) 


where d is the dipole operator defined in Eq. (5 ), and the unit vector € is in the direction of R. The 
matrix element connecting the ground s-state and the lowest, almost degenerate, p-state is calculated 
to be M 0l =(0| d\ l) = -0.395860, and we saw that the energy difference Ao, between these two states 
is 3.996750 x 1 O' 4 Ry, so R 2 is sufficiently small for the perturbation calculation to be valid only 
when /? 2 »1000 or /?»32.This means that for small values of R there are effectively degenerate 
states of opposite parity, coupled by the dipole operator, and hence giving rise to an energy shift that 
is linear in the perturbation, rather than quadratic. (This situation occurs most famously in the 2S 1/2 
and 2P 1/2 states of hydrogen, whose degeneracy is broken slightly by quantum electrodynamic 
effects.) 

A simple way of constructing a variationally correct approximation to the energy shift of the 
ground state (the Stark effect) involves only the two lowest opposite parity states. (We have already 
seen that these produce almost all the asymptotic polarizability.) We write the trial function as 

|T 2 )=C 0 |0> + C,|1>, (it) 


where 0 and 1 refer to the lowest s and p states, respectively. Then the variational expression for the 
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total energy has the form 


£■(«) = 





CiE^ClE^lC^C^MJR 2 


C 0 + Cf 


( 17 ) 


By varying on the two expansion coefficients in the usual way we obtain the eigenvalue condition 


e(R) 


^01 

R 2 


R 2 


^) + A 0 , 


= 0, 


( 18 ) 


which gives rise to a quadratic equation whose solution is 


€(/?) = 



N 


1+4 


M. 


oi 


R 2 A, 


oi / 


-1 


1.998375xl0~ 4 


, 3.924002x1 0 6 , 
1+ -1 


> 

R 4 


( 19 ) 


where we have defined the Stark energy shift of the ground state as e(R) =E Q -E(R ) and have 
inserted the values previously used for M 01 and A 0] . In the large-R limit binomial expansion of the 
square root gives the usual polarizability expression, e(R) = M 0 2 j/(A 0I /? 4 ) , while for small R we 
obtain e(R)=M 0l /R 2 . To display quantitatively the properties of this adiabatic potential, we plot 
a(R) = R*e(R) in Fig. 1. This “field-dependent polarizability’’ goes over to the usual weak-field 
value of a for large R, including only the lowest p-state as discussed above. 

To improve this two-state result in a variationally correct way, we now show two successive 
corrections which introduce the effects of the perturbation of the two lowest-lying states. First, we 
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add to the trial function of Eq. (16) an additional term that takes account of the excited p-states: 



( 20 ) 


where the new term is the perturbation of the ground state, but omitting the lowest p-state. Then we 
can extend this process one more step by adding the perturbation of the lowest p-state: 


5*0 


y(/?)|j)(j|VTO|i) 

E x -E- 


( 21 ) 


To clarify this method, we define in advance all the quantities that will be needed, some of which 
have already been used and give their numerical values. 


M 01 =(l|d|0) = -0.395860, A 01 =3.996750xl0' 4 , a , = £ 


(o|J|^){p|d|o) = 


p* i 


E r E o 


3.207642, 


VE =2.030078, P -E (°l^>(p| d l 0 ) =7.913258. 


s + 0 




hi (E- p ~ E o> Z 


(22) 


P =9 136986 Q = T (QMpXpI^X^I 1 ) _ 6 


s* 0 ( E s~ E l) 


1 (E--e 1 he--e 0 ) 


120350. 


In the evaluation of these sums, 560 intermediate s-states and 728 p-states were retained, with the 
exception of the double sum in Q 3 , where only 120 s-states and 240 p-states were included. 

If we insert the 3- or 4-term trial functions in place of |Y 2 ^ in Eq. (17) and vary on all the 
coefficients, we will obtain a determinantal secular equation for the expectation value of the Stark 
energy shift e(R). For the four-term approximation we find the following form: 
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e{R) 

"o. 

-^1 
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R 2 

R 2 

e(*)+ A 0I 

R 4 
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R 4 


0 

P s e(R) + a s 
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R 4 

R 4 

R 6 

0 


Q 3 


R 4 

R 6 

R 4 


= 0 


(23) 


(This determinant evidently includes the 2x2 and 3x3 subdeterminants appropriate for the other two 
expansion lengths.) In Fig. 2 we have plotted the differences a 3 4 (/?) -oc 2 (/?) between the field- 
dependent polarizability in each of these two expanded approximations and that of the simple 2-term 
value plotted in Fig. 1. As expected from the construction of the 3- and 4-term trial functions, this 
difference goes asymptotically to cc s , the contribution to the polarizability of the ground state coming 
from all intermediate p-states except the lowest one. For other values of R the 4-term result is greater 
than the 3-term result, since they are both variationally correct. 

These curves are the main results of the present work, and they give an effective, R- 
dependent polarizability (or Stark shift) that is a lower bound to the true value. Strictly speaking, 
since this is only a static energy shift it is not directly applicable to the experimental problem of the 
Rydberg states of the HD molecule itself, and it is likely that non-adiabatic corrections and higher- 
multipole contributions will be very large. (We intend to investigate such effects in the future.) But 
it is still interesting to see the effect of the best potential calculated here on a set of Rydberg levels. 
To do this, we simply computed the expectation value 



-13- 


&E NL = -f'¥ flL (R)e l (R)'¥ NL (R)diR, (24) 

using hydrogen wave functions. The results (in MHz) are shown in Table IV; they are, as expected, 
quite large. 
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APPENDIX 

In this work we found it necessary to take account of the near-degeneracy of the two lowest 
states of opposite parity, which contribute the majority of the dipole polarizability in the weak-field 
regime, and are significantly involved in the Stark effect as the electric field increases. One may 
reasonably ask whether higher rotational levels, also of alternating parity, may be important as well. 

The first answer is that since only L=1 states contribute to the perturbation of the L=0 ground 
state, there is no effect on the asymptotic polarizability. But as the field increases and higher-order 
couplings begin to appear these higher-L states should be important. In fact, we intend to extend the 
present research to these higher angular momentum states, but for the moment we describe an 
approximate calculation that shows semiquantitatively that such states are not insignificant. 



-14- 


Let us introduce a rigid rotor model from which we can calculate energy differences and 
dipole matrix elements between the low-lying rotational states. A good approximation is to assume 
that the electron spends about half the time in a hydrogenic state centered on each of the two nuclei. 
Then the rotational energy of each state can be approximated as 



L(L + \) 
2 


( 1 

i 1 

i 

1 1 \ 

1 4 " 




Vp) 


\rn 


1.996745xl0' 4 L(L + l)/ty, 


where fJL DP = 


M D/> (M pD + l) 

m p +m d + 1 


[1224.3435,1224.0104], 


(25) 


and where we have used the data from Table I to evaluate the expectation value and averaged over 
the two configurations. For the energy difference used previously this formula 
gives A 01 =3.993491 x 10“ 4 /?y, in almost exact agreement with the accurate result of Moss that we 
have used in our numerical calculations. We will continue to use this expression for the higher L 
states, although we are aware that there are small centrifugal stretching corrections that should be 
included. 

To calculate the dipole operator in this model, we assume again that there are two terms, and 
that in each one there is a neutral atom (either D or H) and an ion, separated by a distance (^ 12 V The 
contribution from each of these configurations is given by the vector from the center of mass to the 
charged particle. In this way we find the following expression for the dipole operator: 

M -M 

d = p — (r . , \cos6 = 0.684758 cos0 , 

M D +M P + 1 x 12/ 


( 26 ) 
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where we have used the value of ^r [2 ^from Table I, averaged over the lowest s and p states, and 0 
is the angle between the applied electric field and the intemuclear axis. (Incidentally, this formula 
shows explicitly why there is no dipole coupling between the low-lying rotational states for 
homonuclear ions.) Then the dipole matrix element between two neighboring rotational states in this 
model is 

M LLt ,=M L+ i L =°.6 847 58 fdQr^)cos0r° tl (Q) =0.342379 1+ . (27) 

J \ 4L 2 +8L+3 

For the matrix element connecting the two lowest rotational states this formula gives the value 
M m =0.395 345, in excellent agreement with our numerical result. The good agreement both of the 
energy spacing and the dipole matrix element, at least for the first two levels, supports the use of this 
simple model to calculate the Stark effect. 

To evaluate the field-dependent polarizability of HD + with this model, we extend the 
variational trial function of Eq. (16) to include more rotational states. For concreteness we will 
include a total of five. The secular determinant corresponding to Eq. (18) has elements of the 
following general form: 

M 

V[ £ W + V,]V-iHW S V-i]' ls 'V*5. (28) 

K 

By setting the secular determinant to zero we obtain improved values of the energy shift and thus the 
field-dependent polarizability that are almost rigorously variational. In Fig. 3 the effect of including 
first the L=2 state and then the L=3 state is shown. As expected, the polarizability is noticeably 
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increased for intermediate values of R, while the asymptotic polarizability is not affected at all. 
Adding one more term does not make a significant change. 

We expect to find very similar, but more rigorous, results when we carry out a full variational 
calculation including angular momentum terms with L>1 in the trial function. 
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lowest P-state energies and expectation values of various operators for HD + 
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K^&il 

3.143 535 

3.157 374 

HSBjH 

5.486 286 

5.490 379 

■SSI 

5.480 773 

5.484 856 


1.280 997 

1.283 381 

Q b 

1.505 729 

1.509 289 

msm 

- 5.980 024 


mssm 

- 5.982 278 


19 1 

5.385 369 



r.+ r, M D -M r 

r~— +- 


• 1 ’ r -) n pi 

« i 2 . — p — D_f n and z= 


r r. 


12 


2(M p +M D ) 


12 


<2 =— ^r ,2 + r 2 - 3z is the permanent quadrupole moment. 
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Table II. The convergence of the lowest p-state energy and oc^. as a function of Q p , using the 
best-converged form of the ground-state wave function ( Q s =13). 


IE9 

2N a 

E b 


5 

112 

-1.194 238 824 84 

3.089 392 245 14 

6 

168 

-1.194 801 771 34 

3.133 695 375 18 

n 

240 

-1.195 030 443 09 

3.162 882 250 62 

8 

330 

-1.195 230 047 56 

3.184 864 200 87 

9 

440 

-1.195 322 697 52 

3.196 669 457 93 

10 

572 

-1.195 364 513 31 

3.203 457 222 84 

11 

728 

-1.195 371 601 84 

3.207 641 745 15 


a 2N represents the number of terms retained in the expansion. 
b E is the energy of the lowest L=1 state. 
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Table III. The convergence of a p and the ground-state energy as functions of Q s , using the 
best-converged form of the lowest p-state wave function (Q p = 11). 


0, 

N a 

E b 

a p 

7 

120 

-1.195 786 258 30 

2.023 780 342 07 

8 

165 

-1.195 791 941 57 

2.027 828 656 13 

9 

220 

-1.195 794 053 97 

2.028 761 648 12 

10 

286 

-1.195 794 976 73 

2.029 306 062 83 

11 

364 

-1.195 795 415 58 

2.029 653 261 48 

12 

455 

-1.195 795 633 82 

2.029 917 067 13 

13 

560 

-1.195 795 888 91 

2.030 078 303 37 


a N represents the number of terms retained in the expansion 


b E is the energy of the lowest S state. 
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Table IV. Effective potential V HL obtained by averaging e(R) over hydrogen wave functions, and 
the effective polarizability defined by a NL (l/R 4 ')= V NL . 
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FIGURE CAPTIONS 

FIG. I. Field-dependent polarizability a,(/?)of HD + from Eq. (19). The result includes only the two 
lowest states of opposite parity. The dotted line represents the total asymptotic polarizability, 
a = 395.289a 0 3 . 

FIG. 2. Polarizability differences cc 3 (R) - cc 2 (R) (lower curve) anda 4 (/?) - a 2 (R) (upper curve). These 
are solutions of the secular Eq. (23) and its subdeterminants. The dotted line represents the 
asymptotic contribution to the polarizability, a s = 3.2076a 0 3 . 

FIG. 3. Polarizability differences from a 2 (/?)due to the inclusion of higher rotational states in the 
model described in the Appendix. The lower curve includes L=2, and the upper curve also includes 
L=3, and they result from solving the secular determinant with elements given in Eq. (28). Inclusion 
of L=4 has an insignificant effect. 











